rm(list= ls())

library(tidyverse)
library(fixest)
library(fastDummies)
library(lfe)

## Define rank_fun

rank_fun <- function (x) {
  su = sort(unique(x))
  for (i in 1:length(su)) x[x == su[i]] = i
  return(x)
}

## Load Data and Define Outcomes

# Bayern:
# totright2_nofdp = CDU/CSU + AfD 
# totleft2 = SPD and Greens (note that Left is part of other for BY)

bay <- readRDS('data/bavaria_elections_clean.rds') %>%
  mutate(period = rank_fun(year))

## Nr of pre and post treatment periods

id_treat <- bay %>% 
  filter(post == 1) %>% 
  filter(period == min(period)) %>% 
  pull(period) %>% 
  unique()

## Create Period*Treated Indicator Variables 

bay <- fastDummies::dummy_cols(bay, select_columns = 'period') %>% 
  mutate_at(vars(contains('period_')), ~ .*treated) 

## Generate Model Formulas 

outcomevars <- c('totright2_nofdp',
                 'totleft2')
periodvars <- names(bay)[str_detect(names(bay), 'period_')]

## Exclude last pre-treatment period 

periodvars <- periodvars[-which(periodvars == 'period_3')]

results <- lapply(outcomevars, function(dv){
  
  model_formula <- as.formula(paste0(dv, ' ~', paste0(periodvars, collapse = '+'), '| period + AGS |0 | AA.district.id'))
  
  model_est <- felm(model_formula, data = bay) %>%
    broom::tidy(conf.int = T) %>% 
    filter(str_detect(term, 'period')) %>% 
    mutate(outcome = paste(dv))
  
}) %>% reduce(rbind) 

## Rename the period variables

results <- results %>% 
  mutate(period = parse_number(term) - (as.numeric(id_treat) - 1))

## Renaming

results <- results %>% 
  mutate(outcome = recode(outcome, 
                          `totleft2` = 'Total left vote',
                          `totright2_nofdp` = 'Total right vote',
                          `turnout` = 'Turnout')) %>% 
  filter(!str_detect(outcome, '1'))

## Add zero 

zero_period <- data.frame(term = rep('period_2', 2),
                          estimate = rep(0, 2),
                          conf.low = rep(0, 2),
                          conf.high =rep(0, 2),
                          outcome = unique(results$outcome),
                          period = rep(0, 2))

results <- bind_rows(results, zero_period) %>%
  mutate(period = as.factor(period))

## Plot this

myLoc <- 
  (which(levels(factor(results$period)) == "1") +
     which(levels(factor(results$period)) == "0")) / 
  2

p1 <- ggplot(results, aes(x = factor(period), y = estimate*100, 
                          ymin = conf.low*100, ymax = conf.high*100)) +
  geom_hline(yintercept = 0, linetype = 'dotted', color = 'grey40') +
  geom_vline(xintercept = myLoc, linetype = 'dotted', color = 'grey40') +
  geom_errorbar(width = 0) +
  geom_point(fill = 'white', shape = 21) +
  theme_bw() + 
  ylab('Coefficient Estimate\n(percentage points)') + 
  xlab('') + 
  facet_wrap(~ outcome) + 
  scale_color_grey(name = '') +
  scale_x_discrete(breaks = as.factor(seq(-2, 1, 1)),
                   labels = c('2003', '2008', '2013', '2018')) +
  ylim(c(-8, 6))

p1



